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Abstract 

During CO2 sequestration into a deep saline aquifer of finite vertical extent, CO2 will 
tend to accumulate in structural highs such as offered by an anticline. Over times of tens to 
thousands of years, some of the CO2 will dissolve into the underlying groundwater to produce 
a region of relatively dense, saturated water directly below the plume of CO2. Continued 
dissolution then requires the supply of unsaturated aquifer water. In an aquifer of finite 
vertical extent, this may be provided by a background hydrological flow, or a laterally- 
spreading buoyancy-driven flow caused by the greater density of the CO2 saturated water 
relative to the original aquifer water. 

We investigate the long time steady-state dissolution in the presence of a background 
hydrological flow. In steady-state, the distribution of CO2 in the groundwater upstream 
of the aquifer involves a balance between three competing effects: (i) the buoyancy-driven 
flow of CO2 saturated water; (ii) the diffusion of CO2 from saturated to under-saturated 
water; and (iii) the advection associated with the oncoming background flow. This leads to 
three limiting regimes. In the limit of very slow diffusion, a nearly static intrusion of dense 
fluid may extend a finite distance upstream, balanced by the pressure gradient associated 
with the oncoming background flow. In the limit of fast diffusion relative to the flow, a 
gradient zone may become established in which the along aquifer diffusive flux balances 
the advection associated with the background flow. However, if the buoyancy-driven flow 
speed exceeds the background hydrological flow speed, then a third, intermediate regime 
may become established. In this regime, a convective recirculation develops upstream of 
the anticline involving the vertical diffusion of CO2 from an upstream propagating flow of 
dense CO2 saturated water into the downstream propagating flow of CO2 unsaturated water. 
For each limiting case, we hnd analytical solutions for the distribution of CO2 upstream of 
the anticline, and test our analysis with full numerical simulations. A key result is that, 
although there may be very different controls on the distribution and extent of CO2 bearing 
water upstream of the anticline, in each case the dissolution rate is given by the product of 
the background volume flux and the difference in concentration between the CO2 saturated 
water and the original aquifer water upstream. 


1 Introduction 

Carbon capture and storage in deep saline aquifers has been proposed as a potential means to 
limit carbon emissions into the atmosphere, while enabling the continued supply of energy from 
fossil fuels. Much research has been undertaken to explore the processes which control the storage 
of CO 2 over very long periods, and in particular the integrity of a geological storage facility in 
terms of the possible migration of CO 2 back to the surface [a 13 uni US]. Owing to the buoyancy 
of CO 2 relative to water at depths of 1-2 km, CO 2 tends to migrate along permeable sedimentary 
layers and ultimately ponds in structural highs, for example an anticline, which represents the 
upper part of a fold or other deformation in the geological strata (see Figure [Ta]). Such structural 
traps offer a possible storage site providing there is a competent seal rock above the anticline 
(e.g., IPCC[6]). However, CO 2 is soluble in groundwater, which may accommodate concentrations 
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Long time 



(a) Cartoon of the geological problem with region of interest indicated by dashed box. The large red 
arrows represent the direction of the background hydrological flow, while the smaller black curved arrows 
below the CO2, which is trapped at the top of the anticline, represent the convective mixing of the 
CO2 saturated and undersaturated water. The smaller blue arrow inside the dashed box represents the 
buoyancy driven flow of dense CO2 saturated water flowing upstream into the background hydrological 
flow. 
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(b) Model problem for analysis and simulation. 
Figure 1 : Problem of interest. 


of a few wt% CO2 in solution. This in turn leads to an increase in density of the water. With 
the dependency of water density on CO2 concentration, convectively-driven dissolution may 
develop. Water below the trapped CO2 plume becomes increasingly concentrated in CO2 until 
it becomes convectively unstable and sinks into the underlying permeable rock, to be replaced 
by less dense, unsaturated water (c.f. Hewitt et al. [H Lindeberg and Wessel-Berg d Pau et al. 
O Riaz et al. [TI]). Eventually, the water below the CO2 plume becomes fully saturated and the 
continued dissolution requires a more distal supply of undersaturated groundwater. 

Szulczewski et al. m examined the longer-time dissolution by examining the convective 
exchange flow which can develop in a horizontal aquifer. They established that following the 
initial dissolution and near-saturation of the groundwater directly below the plume of CO2, the 
lateral convective exchange flow leads to the slow horizontal spreading of a zone of CO2 enriched 
groundwater associated with the continued dissolution. Eventually, the dynamics of this zone 
may become controlled by a balance between: (a) the buoyancy-driven shear, as the dense 
groundwater spreads along the base of the aquifer; and (b) the vertical diffusion of CO2 from 
this outward spreading dense fluid to the return flow of under-saturated groundwater higher in 
the aquifer. By itself, such buoyancy-driven shear dispersion leads to a progressively waning 
rate of dissolution, and, owing to the relatively low solubility of CO2 in the ground water, the 
prediction that the plume of CO2 may be trapped in the anticline for a very long time m- 

However, at long times the slow background hydrological flows which transport fluid laterally 
through aquifers will become important in controlling the flux of unsaturated water from far 
upstream, especially as other transport processes wane. It is the purpose of this paper to explore 
the long term influence of a background hydrological flow on the process. In this context, 


2 
















Woods and Espie m established some non-linear bounds on the flux of groundwater that may 
reach an anticline along a weakly tilted aquifer resulting from the interaction of a background 
hydrological flow with a convective exchange flow for intermediate times, during which the 
cross-aquifer diffusive transport of CO2 is small. In the present work, we account for the effects 
of such diffusion and this leads to a more complex problem involving the interaction of the 
background advection, the buoyancy-driven flow and the diffusive transport of CO2 upstream of 
the anticline. We note that in our modelling we assume the background hydrological flows are 
constant in time and that there are no mineralogical reactions of the CO2 with the formation; 
these are simplifications but provide a reference with which the effects of mineral precipitation 
or changes in the background forcing over time could be compared. 

We develop a series of idealised, analytical solutions for the governing equations and then test 
these solutions using a full numerical simulation of the two-dimensional governing equations. In 
modelling the flow in porous rock, we assume that the dynamics are governed by Darcy’s Law, 
which relates to slow viscous flow. This is an appropriate model in the present content of slow 
hydrological flows and the slow buoyancy driven flow of CO2 saturated water (e.g. BearfU Woods 
[H]). We thereby establish that when the buoyancy-driven flow of the dense CO2 saturated water 
is large compared to the background flow speed, which is typical, three different regimes may 
become established: (i) the weak diffusion limit in which there is a nearly static intrusion 
of CO2 saturated water upstream; (ii) an intermediate regime in which there is a balance of 
buoyancy-driven flow and vertical diffusion with the oncoming flow; and (hi) a strong diffusion 
limit, in which there is a balance between the upstream diffusion of CO2 and the downstream 
advective transport of unsaturated water. In each case, the dissolution rate is proportional to 
the groundwater flow, even though the controls on the extent of the CO2 enrichment of the 
groundwater upstream of the anticline may be very different. 

2 Model system 

In our analysis, we consider a two-dimensional flow geometry, shown in Figure This corre- 
sponds to an anticline produced by a fold in the geological strata that extends for a relatively 
long distance in the direction normal to the page compared to the width of the fold. We consider 
a background hydrological flow that supplies fluid from the right-hand side, and an aquifer that 
is of uniform thickness and horizontal. We assume that directly below the plume of CO2 which is 
trapped in the anticline, the water is fully saturated in CO2 as a result of the vertical convective 
dissolution, and that this is carried downstream (to the left in Figure [Ta|). The primary purpose 
of this paper is to examine how the concentration of CO2 varies in the upstream direction, shown 
in Figure In our numerical model, we choose the location of the upstream boundary of the 
flow domain to be upstream of the region containing elevated concentrations of CO2, so that we 
can impose a simple uniform flow of fluid with uniform background CO2 concentration. 

The full model involves Darcy flow with a buoyancy term that is dependent on the dissolved 
CO2 concentration c G [cq, cd], where cq > 0 is the initial CO2 mass fraction of the groundwater 
and CD is the mass fraction of the CO2 saturated groundwater below the trapped plume of CO2 
(see Figure [Tbl). We assume a fluid density p given by 

p = po + ^(c-co)po, (1) 

where the constant po ^ 0 is the initial water density and /3 > 0 is the expansion coefficient of 
dissolved CO2 in groundwater. We work with scaled concentration G [ 0 , 1], given by 

c= {cD - Co)c^ -h Cq. ( 2 ) 

To formulate the governing equations in non-dimensional form, we denote dimensionless 
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variables by the superscript and we introduce: 


kI3{cd - co)pog * 

u = - u , 

M 

p = I3{cd- co)pogHp*, 

t = _ ^ _ t\ 

k13{cd - co)pog 

X = Hx*, 

where u is the Darcy velocity, the constants > 0 and /i > 0 are the permeability and viscosity, 
respectively, g is the gravitational acceleration, p is the pressure field, H is the characteristic 
height of the domain, t is time and x is spatial position. 

We denote our domain of interest by D C M^, with boundary V — dO. and outward unit 
normal vector to the boundary n. The boundary is partitioned as depicted in Figure We 
formulate a time-dependent model, with time interval of interest denoted by / = [0,tw)- We 
are interested in steady solutions, hence will be chosen to be suitably large in numerical 
simulations. In terms of non-dimensional quantities, the continuity equation, the Darcy equation 
and the boundary conditions read: 


<1 

0 

on D X /, 

(7) 

= —'V'^p^ + c^Ck 

on D X /, 

( 8 ) 

II 

on Fc X /, 

(9) 

• n — 0 

oil rcap 1 5 

( 10 ) 

n = —u% 

on Fb X /, 

( 11 ) 


where Sk is the unit vector in the direction in which gravity acts, p'^jj is a prescribed pressure 
and > 0 is the prescribed fluid velocity across the inflow boundary. The condition in 0 


gives a background flow from right-to-left in Figure [Tbl 
The concentration of dissolved CO 2 is modelled by: 

dc^ 1 

— + V*c* • u* -V* • V*c* = 0 on X I, (12) 

dt* Ra ^ \ ) 

c*u* ■ n = u* ■ n on Fc^n x I, (13) 

^V*c* • n = 0 on Fc X I, (14) 

Ra 

{-^V*c + c*u*) ■ n = 0 on (Fcap U Fb) X I, (15) 

Ka 

c'^(x,0) = 0 on D, (16) 


(3) 

(4) 

(5) 

( 6 ) 


where Fc,in is the portion of Fc on which • n < 0, the constant Ra is a Rayleigh number, 

k / 3 ( c £) - CQ)pQgH 


Ra — 


fiD 


(17) 


and D > 0 is the pore-scale diffusivity. The boundary condition in (13) ensures that the advective 


flux of CO 2 at the CO 2 trap boundary {x = 0) has dimensionless concentration unity on the 
inflow parts of the boundary, while it is not prescribed on the outflow parts of the boundary. In 


steady state, equations (12)-(15) require that: 


• n ds = 0 


( 18 ) 


We work from this point onward with the non-dimensional equations, hence we drop the 
superscript in the following. 
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3 Physical discussion 


The above non-dimensionalisation identifies two controlling parameters: ub represents the ratio 
of the background flow speed to the buoyancy-driven flow speed, and Ra represents the buoyancy- 
driven flow speed compared to the effective speed associated with vertical diffusive transport 
across the flow domain. These two parameters may be used to delineate the different flow 
regimes which may develop. We explore this below. 


3.1 Gravity intrusion model 

In the case of weak diffusion, we expect that a nearly static intrusion of the dense CO 2 saturated 
fluid extends upstream into the aquifer, and that this is balanced by the pressure gradient of the 
oncoming hydrostatic flow. There will be a thin diffusive boundary layer between the intrusion 
of CO 2 saturated water and the oncoming flow of groundwater. If the intrusion extends far 
into the aquifer (X 1), then the continuity equation suggests that the background flow will 
be largely parallel to the boundary of the domain. To model this regime, we assume a sharp 
interface in the concentration field at a height h{x) above the lower boundary of the aquifer, 
where 0 < h{x) < 1. This interface delineates the CO 2 saturated intrusion and the overlying 
groundwater. Assuming the pressure in the intrusion is approximately hydrostatic (c.f. Huppert 
and Woods [5], Woods Hi)* then in equilibrium the buoyancy driven pressure gradient in the 
x-direction along the intrusion matches the pressure gradient associated with the background 
flow above the intrusion, which has speed ub/— h). This leads to the balance 

The shape of the intrusion is therefore given by 

h{x) = 1 — ^2ubx^ (20) 


(19) 


and it follows that the extent of the intrusion into the aquifer is Xint = 1I2ub- This implies 
that ii Ub is small, the intrusion extends far upstream into the aquifer, relative to the vertical 
extent of the aquifer, and the assumption that the flow is one-dimensional is valid. 

The dimensionless time-of-travel of the oncoming flow past this intrusion is given by Tint = 
^(1/u) dx = l/(3i^^). For the interface to remain sharp, the time should be small relative to 
the diffusion time, Tdiff = Ra. This requires that Ra ^ l/(3u^) which may be expressed in the 
form: 

u%Ra > i. (21) 


For simplicity, we will henceforth use the condition u^^Ra ^ 1. 

When Ub is large, we expect any i ntru sion will become progressively smaller, with Xint < 
1/ub (for example Figure 9a, Section 4.2) and so we now expect the time-scale to be 

an upper bound on the advection time, ^ X[-at/uB- Since Xint < 1 for large we compare 
this with the diffusion time along, rather than across, the aquifer. This diffusion time scales as 
X?^^Ra and suggests that the line i?a = 1 provides an upper bound on the transition from the 


diffusion to the advection regime. We return to this case in Section 4.2 


3.2 Buoyancy-driven shear dispersion model 

In the case <C 1, diffusion in the vertical direction will be relatively fast, hence the vertical 

gradient in concentration across the aquifer will be small. However, there may be a significant 
gradient in the along-aquifer direction. This can lead to different regimes in which diffusion is 
important and we now establish conditions which determine whether a buoyancy-driven shear 
flow develops or a simple advection-diffusion balance controls the transport. We explore these 
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two limits by starting from the full equations and allowing for variations in the concentration of 
CO 2 in the fluid associated with the diffusive flux. We follow largely the analysis of Szulczewski 
et al. m and Woods d to formulate a one-dimensional asymptotic model for the long-time 
evolution of the vertically averaged concentration field, but now in the presence of a background 
flow. 

We decompose the CO 2 concentration of the groundwater in the form c(x,y,t) = c(x,t) -h 
where c is the average concentration across the depth of the aquifer: 


c = 



( 22 ) 


Under the assumptions that the concentration fluctuations c are small, as expected in the limit 
u^^Ra <C 1, and that the horizontal scale of the flow is much larger than the thickness of 
the aquifer, as expected in the case <C 1, the non-hydrostatic vertical pressure gradient is 
relatively small and the flow is approximately parallel to the boundaries of the flow domain. 
Therefore the pressure may be approximated by: 


P^Po- yc, 


(23) 


where po = P^{x->t) is the pressure at the base of the aquifer. 

We decompose the velocity of the fluid u — v) into a sum of the average across the depth 
of the aquifer u = {u^v) and the fluctuation ii = (fi,i)), where: 


f 


udy. 


Using the approximation for the pressure (23), Darcy’s law implies that 

dc 


dpo . - 


and so 


u = 


dc 

dx 


y- 


(24) 


(25) 


(26) 


Taking the vertical average of the transport equation (12), and combining with the continuity 
equation it may be shown that 


dc _dc ^dc du ^ 1 d^c 
dt~^^dx~^^dx~^dx^ Ra dx‘^ ’ 


(27) 


Subtracting (|27|) from the transport equation, we obtain an equation governing the evolution of 

d‘^c 


the concentration fluctuation: 


dc ^ dc _dc ^dc ^dc 1 
+ -^ -^ -^ = 


dt 


dx dx 


dx 


+ 


dy Ra \ dx‘^ dy 


d^c \ ^dc du 
+ u— + —c. 


dx dx 


(28) 


After long time periods, we expect the dominant balance in equation (28) to be between the 
distortion of the mean concentration due to the shear flow and the cross layer diffusion m- 
This gives rise to the following dominant balance, which can be shown a posteriori: 


1 d^c ^ dc 
Ra dy‘^ ^ dx 


(29) 


Inserting (26) into (29) and integrating leads to the expression 
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dc\ y 


y 


1 


c — Rd ( 7-— I I —— — —— -|- —— 


dx J 
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(30) 














Combining this expression with the expression for u in the depth-averaged transport equa- 

1 d^c 


tion (27) becomes: 


dc 


dc 


Ra d 


dt dx 




Ra dx‘^ 120 dx 


dc^ ^ 


(31) 

oc. In the limit that 

u^Ra <C 1, the vertical gradient of concentration is small, and so in this limit it follows from the 


At long times, (31) admits steady solutions in which c 


dx ^ 

0 as X 


boundary condition in (13) that these solutions also require c 1 at x = 0. To help interpret 


these solutions, it is convenient to re-scale the horizontal coordinate according to 

_ ( Ra Y 
~ \120ub) 


X 


X , 


leading to the relation 


where 


dc 


d^c 


__ _ ^ 

dx ^ dx‘^ dx 


dx ) 


a = 


120 

Ra^u^^ 


(32) 


(33) 


(34) 


We see that for large a diffusion dominates (—c = ac'), and for small a dispersion is dominant 
(—c = c'^). We have not found an analytic solution to ( [33| ), but in the two limits a ^ 1 and 
(a <C 1 there are useful analytical approximations. 

In the limit cr <C 1, the buoyancy-driven dispersion balances the advection. The solution 


to (33) when a = 0 is: 




(35) 


Substitution of this solution into equation (28) and comparison of terms identifies that the 
dominant balance is indeed given by (29), in the limit <C 1 and ub < 1- The ffow 


extends a large distance upstream compared to the thickness of the aquifer and the cross-ffow 
diffusion is fast compared to the time for the background ffow to pass through the region in 


which there is an elevated CO 2 concentration. The solution (35) suggests that the region of 


enhanced concentration advances upstream a non-dimensional distance 

Ra \ 3 
2 \120ub ) 


Xdis = 


(36) 


In the limit a ^ 1 the steady-state is dominated by a balance of advection and diffusion, 
and the solution may be approximated by 


- _ ^-RauBX 


with a characteristic length scale of 


^diff — 


In Cdiff 


RauB ’ 

where Cdiff is the concentration at which we consider it to be negligible. 
Equating Xdiff and Xdis, we find that 

3 


OLp — 


(37) 


(38) 


(39) 


2 In Cdiff 

and this provides an indication of the transition between diffusive and dispersive mechanisms. 
For Cdiff = 0.01, we find — 0.326. Figure shows Xdis and Xdiff as function of Ra for 
two different values of the background ffow, using Cdiff = 0.01. For larger Rayleigh numbers, 
dispersion controls the distance that the CO 2 front extends upstream. 
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Figure 2: The distance the CO 2 front extends upstream in the dispersion limit and in the 
diffusion limit as a function of Ra for two different values of background ffow ub> Distances X 
correspond to multiples of the aquifer height H. 


3.3 Regime differentiation 


Combining the analysis of the gravity intrusion with the model of the buoyancy-driven shear 
dispersion, we infer that for ub < 0(1) three regimes may arise, as shown in Figure]^ Gravity 
intrusion occurs when u^^Ra ^ 1. When u^^Ra <C 1, either a diffusion or dispersion dominated 
ffow results, depending on the parameter a (see equation (34)). 

As Ub approaches unity, the above analysis shows that the transition between the dispersion 
and diffusion regimes and between the dispersion and intrusion regimes converge. In the case 
Ub ^ 0(1), the ffow becomes more restricted in lateral extent upstream of the anticline, and our 


analysis of time-scales for > 1 given at the end of Section 3.1, suggests that with ub > 0(1), 
an upper bound for the case in which the along-aquifer diffusion dominates the intrusion regime 
is Ra = 0(1). In the diffusion dominated regime, from the boundary conditions (13) and ( [M] ), 
we expect that at x = 0, c < 1 and we explore this further in Section below. In Figure we 
illustrate this upper bound with a dotted line, which for convenience we show as Ra = 4.93 so 
that it intersects the point at which the solid and dashed lines converge. 


4 Comparison of analytical and numerical models 

To support the asymptotic analysis, the full problem in Sectionj^has been solved on a domain of 
length L = 100 and height H — 1. We use a mixed finite element method for the Darcy ffow, and 
an upwinded discontinuous Galerkin method for the transport equation. Problems are advanced 
in time until a steady-state is reached. A detailed description of the numerical method and the 
complete computer code used to produce all examples is provided in the supporting material HU- 
The computer code is built on the FEniCS libraries [8]. 

4.1 Weak background flows {ub <^) 

To illustrate the form of the velocity and concentration fields for the three different regimes 
when < 1, we show in Figure]^ the computed concentration field for three different values of 
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Figure 3: Illustration of the different regimes. The dashed line denotes a = 1 (see equation (34)) 
which delineates the dispersive and diffusive regimes for ub < 0(1). The solid line shows 
Ra = equation ([2T|)) and delineates the boundary between the intrusion regime and 

(i) the dispersive regime for ub < 0(1) and (ii) the diffusive regime for ub > 0(1). The dotted 
line Ra = 4.93 denotes an upper bound on the transition between the intrusive and diffusive 
regimes in the case ub > when the intrusion is relatively short, and the along-aquifer diffusive 
transport dominates the cross-aquifer diffusion; we have chosen Ra = 4.93, so that this line 
intersects the point where the dashed and solid lines converge. 
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(a) Gravity intrusion: Ra = 3000 
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(b) Dispersive regime: Ra = 50 



(c) Diffusive regime: Ra = 1.0 
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Figure 4: Concentration contours for = 0.1 and different values of Ra. The colours show the 
concentration locally as defined by the scale in the figure and the contour lines are shown at 
equal intervals of 0.1 from 0.1 to 0.9. The x-axis is longer in (c) than (a) or (b) to display the 
full diffusive regime. 


Ra when ub = 0.1. The values of Ra correspond to points in the gravity intrusion, dispersion 
dominated, and diffusion dominated regimes. 


4.1.1 Concentration profiles 

Figure!^ shows the vertically averaged concentration as a function of position along the aquifer 
for the one-dimensional model as given by ( [33| ) (dashed lines) and the two-dimensional full 
numerical calculations (solid lines) for = 8x 10“^ and values of Ra ranging from 100 to 10000. 


The one-dimensional model (equation (33)) was solved using the FEniCS libraries and the full 
code is included in the supporting material [H]. There is very good agreement between the two 
models when Ra < 1000 and the flow is in either the dispersion or diffusion dominated regime. 


The limiting dispersive (35) and diffusive (37) cases (heavy black lines as shown in legend) are 


also shown in the figure, and coincide with the two-dimensional numerical simulation when in 
the appropriate limit regimes. 

When Ra = 10000, the problem is entering the gravity intrusion regime. In Figurej^the two- 
dimensional numerical solution for this case no longer matches the one-dimensional dispersion 


dominated model in (33) since there are significant fluctuations in concentration across the 
height of the aquifer. In this regime the solution has a narrow vertical region of adjustment in 
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Figure 5: Computed and analytical variation of c with x for different values of Ra when ub = 
8 X 10“^. The solid lines represents the two-dimensional numerical solution and the dashed lines 
represent the one-dimensional solution. 
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0 2 4 6 8 10 


Figure 6 : Concentration as a function of position for the gravity intrusion case with = 0.1 
and Ra = 3000 and the prediction of the intrusion model (equation [20|) for the interface position 
overlaid (dotted line). 


the concentration field from the CO 2 saturated fluid at the base of the aquifer to the unsaturated 
oncoming groundwater at the top of the aquifer, reminiscent of the intrusion model given in 
Section 13.11 

Indeed, in Figure we compare the numerical solution for the concentration in the case Ra = 
3000 and 1/5 = 0.1 with the prediction of the interface height h as predicted by equation (20) 


for the gravity intrusion model. That model treats the adjustment of the concentration from 
the CO 2 saturated fluid to the oncoming groundwater flow, as a sharp interface. There is a 
reasonable match for x < 3.5. However, the diffusive boundary layer which is present in the full 
numerical solution leads to a weak recirculation that is not included in the intrusion model. 


4.1.2 Velocity fluctuation profiles 

In the buoyancy driven dispersion regime, which arises for small a and when u^^Ra <C 1, 
equation (26) predicts that the velocity variation from the mean flow u will vary linearly with 


depth. If u is divided by the depth-averaged concentration gradient, the profiles are predicted 
to pass through —0.5 at y = 0 and 0.5 at y = 1. Figure [7a| shows the scaled velocity profiles 
computed from the two-dimensional model for Ra = 1000 and ub = 8 x 10“^ {a — 0.0527). For 
X < {8/A)X(iis the velocity profiles are linear, whereas close to the stall point {x/ XfUs — 1) the 
profile begins to deviate from the simplified theory. 
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(a) Velocity fluctuations u scaled by concentration gradient for Ra = 1000 
and UB = Sx 10“^ {a = 0.0527). 



(b) Velocity fluctuations u for different values of a at different values of x. 

Figure 7: Scaled velocity profiles at different points along the domain for different values of a 
with ub = 8 X 10 “^. 
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Figure 8: Typical streamlines of the flow for = 0.01 and Ra — 1,10 and 100. 


When in the regime where diffusion is important {u^^Ra <C 1), as a increases and the 
flow transitions from dispersive to diffusive, the one-dimensional analytical model becomes less 
applicable for u. This can be seen in Figurewhere the scaled horizontal velocity fluctuations 
are shown at different distances into the domain for different values of a. For 0.0572 < a < 
0.1442, the scaled horizontal velocity fluctuations at different distances into the domain all lie on 
a straight line between —0.5 and 0.5. As a increases, the numerically computed profiles deviate 
from the model, as the lateral diffusive transport becomes comparable to and then progressively 
more significant than the shear dispersion, but in this diffusion limit these perturbation velocities 
are very small compared to the background hydrological flow. When Ra = 10000 and ub = 8 x 
10“^ {a = 0.0027) the gravity intrusion region is being approached and the vertical concentration 
fluctuations across the domain are no longer small and so the velocity fluctuations are no longer 
governed by equation (26). 

In Figure]^ we present a picture of the typical streamlines for the flow when ub = 0.01 and 
Ra = 1,10 and 100. The streamlines are computed, approximately, by solving equation (33) 
numerically and inserting the result into (26), and then adding on the background hydrological 
flow. The case i?a = 1 corresponds to the diffusion limit, the case Ra = 100 corresponds to the 
dispersion regime, and the case Ra = 10 is an intermediate case. The figure shows that in the 
diffusion limit the oncoming groundwater flow is dominant and the flow remains nearly uniform, 
but as Ra increases, some fluid begins to flow upstream leading to a small circulation in the lower 
part of the aquifer. As Ra increases further, and the flow is controlled by the buoyancy driven 
dispersion, a strong recirculation develops upstream of the anticline, leading to the diversion of 
the oncoming groundwater flow towards the top of the aquifer. 


4.2 Strong background flows {ub ^ 1) 

We now look at the case with a stronger background flow {ub ^ 1). Three examples of full 
numerical solutions are shown in Figure [^corresponding to Ra = 1000, 2 and 0.5 for ub = 1.0. 
In contrast to the case ub < 0(1), only two distinct regimes develop (figure 3) . Now, the 
flow transitions from the gravity intrusion regime to the diffusion dominated adjustment of 
the concentration since, with large ub^ the upstream extent of the buoyancy driven flow is 
insufficient for the buoyancy driven dispersion to develop before the along aquifer diffusion 
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(a) Gravity intrusion regime: Ra = 1000 
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(c) Diffusive regime: Ra = 0.5 
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Figure 9: Concentration fields for different values of Ra for = 1.0. The white contour in (b) 
and (c) is the c = 0.01 contour. The presented domains have been truncated at different lengths 
to best show each regime. Note that the maximum concentration for the colour bar is set equal 
to 0.3, which is greater than the highest concentration for panels (b) and (c). 
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Figure 10: Computed vertically averaged concentration for various values of Ra and at x = 0. 


becomes significant. With large i?a, the solution is similar to the gravity-driven intrusion solution 
(Figure [9a|), while for smaller values of Ra the solution evolves towards the along aquifer diffusion 
solution (Figure!^. 

We have calculated the vertically averaged mean concentration at x = 0 from the numerical 
solutions for the cases when Ra — 10,1,0.1 and 0.01, as shown in Figure [T0| As ub increases, c 
decreases at x = 0 since there is a progressively stronger flow from the upstream region which 
suppresses the upstream buoyancy driven flow of dense, CO 2 saturated fluid from below the 
anticline 

In the diffusive regime, which we expect to apply for small i?a, the vertically averaged 


concentration may be approximated by the diffusive solution of equation (33), given by: 


cix) = 


(40) 


Figure shows the vertically-averaged concentration profiles for the three values of Ra and 
Ub — 1.0. Using the numerically determined value for the mean concentration at x = 0, 
we have compared the vertically averaged concentration with the diffusion solution given by 
equation ( [40| ). When Ra = 0.1, the system is in the diffusion regime and the numerical solution 
for c matches the diffusion profile. When Ra > 1, the simulations move towards the intrusion 
regime and the numerical solutions evolve away from the approximate analytical solution. 


5 Conclusions 

We have explored both analytically and numerically the long-term dissolution of a plume of 
CO 2 trapped in an anticline and driven by a steady background flow of CO 2 unsaturated water 
from upstream. We have focused on the role of diffusion and buoyancy-driven flow in regulating 
the distribution of CO 2 in the aquifer fluid upstream of the anticline. In the case ub < 
where is the dimensionless background hydrological flow, the buoyancy-driven speed of the 
dense CO 2 saturated water exceeds the oncoming flow speed and the CO 2 extends a significant 
distance upstream of the anticline {x ^ H). In this case, we have established that three different 
regimes may develop. With a small diffusive flux across the aquifer {u^^Ra ^ 1, where Ra is 
the Rayleigh number associated with the dense CO 2 laden fluid), a static intrusion of dense 
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Figure 11: Computed and analytical variation of c with x for different values of Ra when ub = 
1.0. The solid line represents the two-dimensional numerical solution and the dashed line repre¬ 
sents the diffusion limit if using the numerically computed c as the boundary condition at x = 0. 


CO 2 saturated fluid develops and extends a distance Xj^UB upstream of the anticline. This is 
balanced by the pressure gradient in the oncoming flow. With larger diffusive fluxes across the 
aquifer {u^^Ra <C 1) we have established that a buoyancy-driven shear dispersion flow regime may 
develop and a convective recirculation develops just upstream of the aquifer, regulated by (i) the 
supply of unsaturated aquifer fluid from upstream; (ii) the buoyancy-driven flow associated with 
the dense CO 2 saturated fluid from downstream; and (iii) the vertical diffusion of CO 2 across 
the aquifer. However, if the diffusive transport is too rapid (cr <C 1) then a simple advection- 
diffusion balance regulates the distribution of the CO 2 in solution in the water upstream of the 
anticline. In the case ub > the CO 2 extends a much smaller distance upstream from the 
anticline, and in this case, either only the advection-diffusion balance or the intrusion regimes 
develop. 

In the context of CO 2 sequestration in deep saline aquifers, this analysis is important as 
it demonstrates the strong effect that a background hydrological flow has on the long-term 
dissolution of CO 2 in a structural trap. We now show that under some typical conditions the 
dynamics may indeed be controlled by a balance between the buoyancy-driven shear dispersion 
and the background hydrological flow. This leads to new estimates of the maximum upstream 
migration of CO 2 rich groundwater. The solubility of CO 2 in groundwater is only a few wt%, so if 
we consider a plume of CO 2 of order 10 m deep, trapped in a structural anticline and connected to 
a laterally extensive aquifer of order 20-30 m deep, then vertical convective dissolution alone will 
only lead to dissolution of order 0.2-0.6 m. Continued dissolution will require the lateral supply 
of undersaturated water from the aquifer and this may be achieved through a combination of 
buoyancy-driven lateral dispersion of the dense CO 2 saturated water from below the CO 2 plume 
and supply of water resulting from a background hydrological flow. For typical conditions, 
with permeability of order 0.1-0.01 Darcy, a density difference of order a few percent between 
the undersaturated and saturated water, and an aquifer diffusivity of order 10“^-10“^^ m^/s, 
the Rayleigh number will be of order 10^-10"^. With a hydrological flow speed of order 10“^- 
10“^ m/s, the dimensionless velocity ub will be of order 0.01-1.0. From Figure]^ we see that it 
is the transport associated with the shear dispersion that balances the steady background flow. 
We estimate that the length scale of the dispersive transport (equation ([3^) will be of order 
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100-400 m. Once the steady flow regime is established, the continued dissolution will occur at 
a rate proportional to the supply of undersaturated water in the hydrological flow, as given in 
non-dimensional form by ub{cd — cq)- For an anticline whose extent in the direction of the flow 
is of order 1000 m, and with a CO 2 plume with initial depth of order 10 m, then in order to 
dissolve, this will require a net flow of groundwater of order 10^ m^, which will require a time 
of order lO^^-lO^^ s corresponding to 10^-10^ years. 
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